package Modules::JunctionRW;
use strict;
use Data::Dumper;
use Carp;
use Bio::DB::Sam;
sub new {
  my ($class,$util) = @_;
  my $self = {
               util => $util,
             };
  bless $self,$class;
  return $self;
}

sub parse_sam {
  my $self = shift;
  my $util = $self->{util};
  my $path = $util->{DATA_PATH};
  my $file = $util->{SAM_FILE };
  my $filename  = $util->{DATA_PATH}.$util->{SAM_FILE};
  my ($fho1,$fho2,$fho3,$junction_reads);

  open ($fho1,'>'.$path.$file.'.spliced_reads')      or die "$path$file.spliced_reads $!";
  open ($fho2,'>'.$path.$file.'.spliced_reads_info') or die "$path$file.spliced_reads_info $!";
  open ($fho3,'>'.$path.$file.'.unspliced_reads')    or die "$path$file.unspliced_reads $!";

  my $bam          = Bio::DB::Bam->open($filename);
  my $header       = $bam->header;
  my $target_count = $header->n_targets;
  my $target_names = $header->target_name;
  
  while (my $align  = $bam->read1) {
    my (@matches,@introns);
    my $start            = $align->pos+1; 
    my $cigar            = $align->cigar_str;
    my $chr              = $target_names->[$align->tid];
    next if $cigar eq '*'; #no alignment
    Carp::croak "unexpected cigar string: $cigar in ".$align->qname."\n" if ($cigar=~/[^\d+|M|N|I|D]/g);
       
    if(@matches = ($cigar=~/([^N]*)(N)([^N]*)/g)) {
      $junction_reads++;  
      my ($locator,$flank_start,$flank_end,$junc_start,$junc_end,$no_mismatches);
      print $fho1 $align->data."\n";
      my $pruned_matches = "@matches";         
      my @matches        = split (/\s+/,$pruned_matches);
      my $match_counter  = 0;
               
      foreach my $match (@matches) {     
        if ($match eq 'N' ) {          
          $flank_start        = $junc_end || $start;
          $locator            = $junc_end || $start;
          my $left            = $match_counter-1; 
          my $index           = $match_counter;
          my $right           = $match_counter+1;
          my $left_flank      = $matches[$left ];
          my $right_flank     = $matches[$right];                
          my $junction        = $left_flank.$match.$right_flank; 
                
          ($junc_start,$junc_end,$no_mismatches) = $self->calc_read_locator(\$locator,$junction);

          $flank_end          = $locator;  

          print $fho2 "flank_start:     ".$flank_start;                 
          print $fho2 "\tjunc_start:    ".$junc_start;
          print $fho2 "\tjunc_end:      ".$junc_end;         
          print $fho2 "\tflank_end:     ".$flank_end;
          print $fho2 "\tno_mismatches: ".$no_mismatches;
          print $fho2 "\n";
            
          $self->calc_junctions($junc_start,$junc_end,$flank_start,$flank_end,$no_mismatches,$align,$chr);          
        }
        $match_counter++;
      } 
    #last if $junction_reads == 10;
    } else {
       print $fho3 $align->data."\n";
    }  
  }
  print "No. of read junctions:        " . $junction_reads . "\n";
  print "No. junctions w/o filtering:  " . scalar keys %{$self->{'junctions'}};
  print "\n";

  close $$fho1,$fho2,$fho3;
  return $self->{'junctions'};
}
sub calc_junctions {
  my ($self,$junc_start,$junc_end,$flank_start,$flank_end,$no_mismatches,$align,$chr) = @_;
  my $util                   = $self->{'util'};
  my $path                   = $util->{'DATA_PATH'};
  my $file                   = $util->{'SAM_FILE' };
  my $MAX_MAPPINGS_PER_READ  = $util->{'MAX_MAPPINGS_PER_READ'};
  my $junction;
  my $id_junction            = $junc_start."_".$junc_end."_".$chr;
  my $read_name              = $align->qname; 
  my $flag                   = $align->flag;
  my $strand                 = $align->strand;
  my $QUAL                   = $align->qual;        # The quality score for the alignment as a whole
  my $tid                    = $align->tid;         # Return the target ID of the alignment
  my $mtid                   = $align->mtid;        # For paired reads, the target ID of the mate's alignemnt
  my $mpos                   = $align->mpos;        # For paired reads, the 0-based leftmost coordinate of the mate's alignment on the reference sequence
  my $cigar                  = $align->cigar_str;
  my $is_paired              = $align->paired;      # 1 if the aligned read is part of a mate/read pair 
  my $is_proper              = $align->proper_pair; # 1 if the aligned read is part of a mate/read pair and both partners mapped to the reference sequence
  my $is_unmapped            = $align->unmapped;    # 1 if the read failed to align
  my $mate_is_unmapped       = $align->munmapped;   # 1 if the read's mate failed to align
  my $non_redundant_reads    = join "_" ,$flank_start,$junc_start,$junc_end,$flank_end; 
     
  push @{$self->{'junctions'}->{$id_junction}{'left_flank'   }},$flank_start;
  push @{$self->{'junctions'}->{$id_junction}{'right_flank'  }},$flank_end;
  push @{$self->{'junctions'}->{$id_junction}{'no_mismatches'}},$no_mismatches; 
  push @{$self->{'junctions'}->{$id_junction}{'read_names'   }},$read_name;
  push @{$self->{'junctions'}->{$id_junction}{'strand'       }},$strand; 

  $self->{'junctions'}->{$id_junction}{'support'            }++;
  $self->{'junctions'}->{$id_junction}{'junc_start'         } = $junc_start;
  $self->{'junctions'}->{$id_junction}{'junc_end'           } = $junc_end;
  $self->{'junctions'}->{$id_junction}{'chr'                } = $chr;                                                                 
  $self->{'junctions'}->{$id_junction}{'paired'             }++ if $is_paired;
  $self->{'junctions'}->{$id_junction}{'proper'             }++ if $is_proper;
  $self->{'junctions'}->{$id_junction}{'unmapped'           }++ if $is_unmapped;
  $self->{'junctions'}->{$id_junction}{'mate_unmapped'      }++ if $mate_is_unmapped;
  $self->{'junctions'}->{$id_junction}{'non_redundant_reads'}{$non_redundant_reads}++; 

  return;
}
sub calc_read_locator {
  my ($self,$locator,$flank) = @_;
  my ($no_mismatches,$junc_start,$junc_end) = (0,$$locator,$$locator);
  $$locator--;
  while ($flank=~/(\d+)([M|N|I|D])/g) {
    my ($size, $type) = ($1, $2);
    if ($type eq 'M') { $$locator   += $size }
    if ($type eq 'I') { $$locator    = $$locator - $size +1;
                        $no_mismatches++;   }
    if ($type eq 'D') { $$locator    = $$locator + $size ; 
                        
                        $no_mismatches++;   }
    if ($type eq 'N') { 
                        $junc_start = $$locator;
                        $$locator   += $size;
                        $junc_end   = $$locator +1;                    
    } 
  }
  return ($junc_start,$junc_end,$no_mismatches);
}

1;